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Abstract. We study the quantum mechanics of a billiard (Robnik 1983) in the regime 
of mixed-type classical phase space (the shape parameter A = 0.15) at very high-lying 
eigenstates, starting at about 1.000.000th eigenstate and including the consecutive 
587654 eigenstates. By calculating the normalized Poincare Husimi functions of the 
eigenstates and comparing them with the classical phase space structure, we introduce 
the overlap criterion which enables us to separate with great accuracy and reliability 
the regular and chaotic eigenstates, and the corresponding energies. The chaotic 
eigenstates appear all to be dynamically localized, meaning that they do not occupy 
unformly the entire available chaotic classical phase space component, but are localized 
on a proper subset. We find with unprecedented precision and statistical significance 
that the level spacing distribution of the regular levels obeys the Poisson statistics, 
and the chaotic ones obey the Brody statistics, as anticipated in a recent paper by 
Batistic and Robnik (2010), where the entire spectrum was found to obey the BRB 
statistics. There are no effects of dynamical tunneling in this regime, due to the high 
energies, as they decay exponentially with the inverse effective Planck constant which 
is proportional to the square root of the energy. 
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1. Introduction 

In quantum chaos [H EJ [3] of general (generic) time-independent (autonomous) Hamilton 
systems in the strict semiclassical limit we can conceptually separate regular and chaotic 
eigenstates. This picture goes back to the work by Percival in 1973 pE] and Berry 
and Robnik in 1984 [5]. One of the main results in quantum chaos is the fact that 
in classically fully chaotic (ergodic, autonomous Hamilton) systems with the purely 
discrete spectrum the fluctuations of the energy spectrum around its mean behaviour 
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obey the statistical laws described by the Gaussian Random Matrix Theory (RMT) 
[6j [7], provided that we are in the sufficiently deep semiclassical limit. The latter 
condition means that all relevant classical transport times, like the typical ergodic time, 
or diffusion time, are smaller than the so-called Heisenberg time, or break time, given 
by tn = 2nh/AE, where h = 2ttH is the Planck constant and AE is the mean energy 
level spacing, such that the mean energy level density is p(E) = 1/ AE. This statement 
is known as the Bohigas - Giannoni - Schmit (BGS) conjecture and goes back to their 
pioneering paper in 1984 [8j, although some preliminary ideas were published in [9]. 
Since AE oc h d , where d is the number of degrees of freedom (= the dimension of the 
configuration space), we see that for sufficiently small fi the stated condition will always 
be satisfied. Alternatively, fixing the H, we can go to high energies such that the classical 
transport times become smaller than tn- The role of the antiunitary symmetries that 
classify the statistics in terms of GOE, GUE or GSE (ensembles of RMT) has been 
explained in [10], see also [TT] and [TJ [2j [3j |6]. The theoretical foundation for the 
BGS conjecture has been initiated first by Berry [12], using the Gutzwiller periodic 
orbit theory (trace formula) [13] (for an excellent exposition see |lj) and later further 
developed by Richter and Sieber [TJ], arriving finally in the almost-final proof proposed 
by the group of F. Haake [US OH QU [18] . 

On the other hand, if the system is classically integrable, Poisson statistics applies, 
as is well known and goes back to the work by Berry and Tabor in 1977 (see [H [2, |3] 
and the references therein, and for the recent advances |19j). 

In the mixed-type regime, where classical regular regions coexist in the classical 
phase space with the chaotic regions, being a typical KAM-scenario which is the generic 
situation, the so-called Principle of Uniform Semiclassical Condensation (of the Wigner 
functions of the eigenstates; PUSC) applies, based on the ideas by Berry [20], and further 
extended by Robnik [3]. If the stated semiclassical condition is satisfied, the chaotic 
eigenstates are uniformly extended, and consequently the Berry- Robnik statistics [21] 
is observed - see also [3]. If the semiclassical condition stated above requiring that tn 
is larger than all classical transport times is not satisfied, the chaotic eigenstates will 
not be extended but localized and the Berry-Robnik statistics must be generalized as 
explained in [221 US El [23 [25] and in this paper. 

The relevant papers dealing with the mixed-type regime after the work [5] are [21] 
- [31] and the most recent advance was published in |24j. If the couplings between the 
regular eigenstates and chaotic eigenstates become important, due to the dynamical 
tunneling, we can use the ensembles of random matrices that capture these effects 
[321 [24] . As the tunneling strengths typically decrease exponentially with the inverse 
effective Planck constant, they rapidly disappear with increasing energy, or by decreasing 
the value of the Planck constant. In this work we shall deal only with high- lying 
eigenstates, and therefore we can neglect the effects of tunneling. 

However, quite generally, if the semiclassical condition is not satisfied, such that tn 
is no longer larger than the relevant classical transport time, like e.g. the diffusion time in 
fully chaotic but slowly ergodic systems, we find the so-called dynamical localization, 
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or Chirikov localization. Dynamical localization was discovered in time dependent 
systems [33]. It was intensely studied since then in particular by Chirikov, Casati, 
Izrailev, Shepelyanski and Guarneri, in the case of the kicked rotator as reviewed in 
[31]. See also the references [35]- [38], and the most recent work [39]. For a general 
overview of the time dependent Floquet systems see also [HE]- It has been observed 
that in parallel with the localization of the eigenstates one observes the fractional 
power law level repulsion (of the quasienergies) even in fully chaotic regime (of the 
finite dimensional kicked rotator), and it is believed that this picture applies also to 
time independent (autonomous) Hamilton systems and their eigenstates [39]. (See the 
excellent review of localization in time independent billiards by Prosen in [10].) Indeed, 
this has been analyzed with unprecedented precision and statistical significance recently 
by Batistic and Robnik |24J in case of mixed-type systems, and the present work is being 
extended in the analysis of separated regular and chaotic eigenstates. An early attempt 
of separation of eigenstates in the billiard system has been published in [U], using a 
different approach at much lower energies and with much smaller statistical significance. 

In this paper we introduce a criterion for classifying eigenstates as regular and 
chaotic, and moreover, we show that the regular levels obey the Poisson statistics, 
whilst the chaotic dynamically localized eigenenergies obey exceedingly well the Brody 
distribution |42] , with the Brody parameter values (3 within the interval [0,1], where 
(3 = yields the Poisson distribution in case of the strongest localization, and = 1 gives 
the Wigner surmise (2D GOE, as an excellent approximation of the infinite dimensional 
GOE), which describes the extended chaotic eigenstates. It turns out that the Brody 
distribution introduced in |42j, see also [13], fits the empirical data much better than 
e.g. the distribution function proposed by F. Izrailev (see [35J [3l] and the references 
therein) . 

It is well known that Brody distribution so far has no theoretical foundation, but 
our empirical results show that we have to consider it seriously in dynamically localized 
chaotic eigenstates, thereby being motivated for seeking its physical foundation, and an 
analogous result was obtained in the recent work of Manos and Robnik (2013) [39] in the 
analysis of the quantum kicked rotator, where the object of study are the eigenstates 
of the Floquet operator and the statistical properties of the spectrum of eigenphases 
(quasienergies) in classically fully chaotic regime. 

In the Hamilton systems with classically mixed-type dynamics, which is the generic 
case, we have classically regular quasi-periodic motion on <i-dim invariant tori (d is the 
number of freedoms) for some initial conditions (with the fractional Liouville volume 
pi) and chaotic motion for the complementary initial conditions (with the fractional 
Liouville volume P2 = 1 — pi). The chaotic set might be further decomposed into several 
chaotic regions (invariant components) in case d = 2, whilst for d > 2 it is strictly 
speaking always just one chaotic set due to the Arnold diffusion on the Arnold web, 
which pervades the entire phase space. In sufficiently deep semiclassical limit the Berry- 
Robnik picture [5] is established, based on the statistically independent superposition 
of the regular (Poissonian) and chaotic level sequences, based on PUSC, as explained 
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above. 

In this paper we shall consider only the case of just one chaotic component, although 
the results can be easily generalized for more than one chaotic component. This picture 
gives an excellent approximation for the statistics of spectral fluctuations of the mixed- 
type systems, if the largest chaotic component is much larger than the next largest one, 
which typically indeed is the case e.g. in 2D billiards. 

The energy spectrum of the mixed-type system with one regular and one chaotic 
component can be described in the Berry-Robnik (BR) regime of sufficiently small 
effective Planck constant H e ff by the following formula for the gap probability E(S), 

E(S) = E r ( Pl S)E c (p 2 S) (1) 

and the level spacing distribution P(S) (see e.g. [3]) is of course always given as the 
second derivative of the gap probability, namely P(S) = d 2 E(S)/dS 2 , so that we have 

P(S) = ^E M S)EMS) = ^ + 2§ § + E**. (2) 

This factorization formula (DQ) is a direct consequence of the statistical independence, 
justified by PUSC. Here by E r (S) = exp(— S) we denote the gap probability for the 
Poissonian sequence with the mean level density one. By E C (S) we denote the gap 
probability for the chaotic level sequence with the mean level density (and spacing) one. 
Note that the classical parameter pi and its complement p 2 = 1 — pi enter the expression 
as weights in the arguments of the gap probabilities. 

Using the Bohigas-Giannoni-Schmit conjecture we conclude that in the sufficiently 
deep semiclassical limit E C (S) is given by the RMT, and can be well approximated by 
the Wigner surmise 

MS) = ^exp(-^V F w (ST) = l-W w {Sr) = e xp(-^\ (3) 

such that E C (S) is equal to 

^(S) = 1 - erf (^)= eric (^), (4) 

where erf(x) = ^= J* e~ u2 ' du is the error integral and erfc(x) its complement, i.e. 
erfc(a;) = 1 — erf(x). In the equation (J3J) Ww(S) denotes the cumulative Wigner level 
spacing distribution Ww(S) = J Pw( x ) dx and Fw its complement. The explicit 
Berry-Robnik level spacing distribution (in the special case of one regular and one 
chaotic component) follows immediately, 



(5) 



The correctness of this distribution function in the BR regime (sufficiently small h e ff) 
is by now very well established in highly accurate numerical calculations for all E(k, L) 
probabilities, not only the gap probability [2Tj . 
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In the present work the above basic BR formula ([T]) is generalized as in [21] to 
capture the dynamical localization effects, reponsible for the deviation from the BR 
regime. 

At not sufficiently small h e ff (e.g. in billiards this means at low energies) the chaotic 
eigenstates (their Wigner functions in the phase space) are not uniformly extended over 
the entire classically allowed chaotic component, but are dynamically localized. Thus we 
see the transition from GOE in case of extended chaotic states to the Poissonian statistics 
in case of strong localization. The level spacing distribution in such a transition regime 
of localized chaotic eigenstates can be described by the Brody distribution with the only 
one family parameter /3, 

P B (S) = d^exp (-C 2 ^ +1 ) , W B {S) = 1 - exp (-C 2 S P+1 ) , (6) 

where the two parameters C\ and C% are determined by the two normalizations 
< 1 >=< S >= 1, and are given by 

c 1= ^+i)c 2 , ^=( r (f$i)y +i ( ? ) 

with T(x) being the Gamma function. As mentioned before, if we have extended chaotic 
states = 1 and RMT ([3]) applies, whilst in the strongly localized regime /3 = and 
we have Poissonian statistics. Again, by Wb(S) we denote the cumulative Brody level 
spacing distribution, Wb(S) = Pb(x) dx. The corresponding gap probability is 

where Q(a,x) is the incomplete Gamma function 

POO 

Q(a,x) = / e^e-'dt. (9) 

J X 

By choosing E C (S) in equation (JTJ as given in (JBJ) we are able to describe the localization 
effects on the chaotic component. Such approach has been already proposed in the paper 
by Prosen and Robnik in 1994 [22| 123] and the resulting level spacing distribution, 
emerging from this assumption, was called Berry- Robnik-Brody (BRB). It has two 
parameters, the classical parameter p\ and the quantum parameter /3. It will turn 
out that this description is indeed excellent, and has been verified to high accuracy by 
Batistic and Robnik |24j . 

The applicability of the Brody distribution in this context is theoretically still not 
well understood, but we shall see that the theory describes very well the empirical 
data from the highly accurate energy spectra of billiards at energies around and below 
the deep semiclassical (Berry-Robnik) regime. Therefore, the resulting theory is of 
semiempirical nature, but seems to be quasi-universal in the sense that below the 
Berry-Robnik regime we indeed find spectral fluctuations, in particular the level spacing 
distribution, which are well described by our theory on the very finest scale of level 
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spacings. We have also tried to use other one-parametric level spacing distributions 
instead of Brody, in particular those proposed and studied in Izrailev's papers [35 | l38 t 134"] 
(see also [33]), but must definitely conclude that the Brody distribution is quite special, 
it gives by far the best agreement between the theory and real spectra. Izrailev's 
intermediate level spacing distribution intended to capture the localization effects 
manifested in the the quantal spectra is given by 



where the constants A and B are determined by the normalizations < 1 >=< S >= 1. 

As we shall show the dynamical localization effects can persist up to very high-lying 
eigenstates, even up to one million, whilst - as mentioned before - the tunneling effects 
occur usually only at very low-lying eigenstates, due to the exponential dependence on 
the reciprocal effective Planck constant, oc exp(— const ./ h e f /) , and thus can be neglected 
in our case. 

The paper is structured as follows: In section [2] the billiard system is defined 
as introduced by Robnik (35j HE] with shape parameter A = 0.15, and we describe 
the Poincare Husimi functions, in section [3] we introduce the method of classifying 
and separating regular and chaotic eigenstates in terms of normalized Poincare Husimi 
functions (which are Gaussian-smoothed Wigner functions), in section H] we present the 
results, and in section [5] we conclude and discuss the results. 

2. Introducing the model system and the definition of the problem 

As an interesting and frequently studied model system we have chosen the billiard 
introduced by Robnik in references [15]- [16], whose boundary is defined by the quadratic 
conformal map of the unit circle \z\ = 1 of the ^-complex plane onto the w-complex plane 
(which is the physical plane) as follows 



The choice of this quadratic map has two reasons: (i) it allows for an elegant method [46J 
to solve the Helmholtz equation in the w-plane by transforming back to the z-plane, and 
(ii) it is the simplest one with nontrivial classical dynamics [45J. The shape parameter 
A goes from (circle; integrability) to 1/2 (cardioid billiard; ergodicity: full chaos [4"T]). 
For < A < 1/4 the billiard boundary is convex and thus we observe the existence 
of Lazutkin's caustics [HI [49] near the boundary in the configuration space, which are 
the projections of the KAM invariant tori (in the phase space). When increasing the 
value of A from we have at A = 1/4 for the first time a point of zero curvature 
located at z = —1, and thus w = — 1 + A. By a theorem due to J. Mather this is 
a sufficient condition for the destruction of Lazutkin's caustics and of the underlying 
invariant tori near the boundary, which in turn is a necessary condition for the ergodicity 
of the classical billiard dynamics. However, at A > 1/4 the dynamics is not yet ergodic, 




(10) 



w = z + Xz 2 . 



(11) 
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fully chaotic, as there are still some tiny KAM islands of stability |50j, not so easy to 
detect numerically. Nevertheless, for A = 1/2, the ergodicity was proven rigorously by 
Markarian |47_j. 

We are interested in the mixed-type case, within the interval < A < 1/4, 
where the chaotic regions coexist with the regular regions in the phase space. In 
particular, we have chosen A = 0.15, which is one of the most frequently studied cases 
in [261 [271 [221 [231 [251 EH [521 [H] . The fractional phase space volume p 1 of the classically 
regular part of the phase space (not to be confused with the area on the Poincare surface 
of section!) is equal to 0.175, as has been carefully studied in [511 [21] . In a previous 
work [221 [23] p\ was estimated numerically as p\ = 0.36, which is due to the technical 
difficulties in distinguishing the regular regions and slowly diffusing chaotic regions due 
to the sticky objects in the classical phase space. These difficulties were overcome in [51] 
and |24j, using the new methods based on the ideas and approach in [53]. The estimate 
Pi = 0.175 is now believed to be accurate within at least one percent relative error. 

The classical mechanics of 2D billiards is studied in the Poincare-Birkhoff 
coordinates (s,p), where s is the arclength parameter going from to £, the perimeter 
of the billiard domain, in our case counted anticlockwise from the point w = 1 + A, 
and p is simply the sine of the reflection angle, which goes from -1 to +1. The bounce 
map is defined by the free motion between the collision points on the boundary, obeying 
the specular reflection law upon each collision. Thus, the complete information on 
the classical dynamics is contained in the structure of the bounce map on the cylinder 
(s,p). When analyzing the quantum mechanics of this system, we would like to find an 
analogous two dimensional space which also contains the complete information about 
the quantum mechanics, namely about the eigenfunctions. This is the space of the so- 
called Poincare Husimi functions (see [51] and the references therein) that we introduce 
below. 

The quantum mechanics of the billiard system comprises the study of the solution 
of the Helmholtz equation for the billiard domain B, 



with the Dirichlet boundary condition ip = on the boundary OB. We have used a 
number of methods, like in [21], to calculate the eigenenergies Ej = fc?, where j is the 
counting index j = 1, 2, 3, . . ., and the associated eigenfunctions are denoted by ipj(r). 

Introducing the important quantity u(s), the normal derivative of the eigenfunction 
ip on the boundary, from here onwards called boundary function, 



where n is the unit outward vector normal to the boundary at position s, we can show 
[55] that the eigenvalue problem f[T2|) is equivalent to the following integral equation 




Aip + k 2 ip = 0, 



(12) 



u(s) = n • V r ip (r(s)) , 



(13) 





(14) 
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Here r(t) is the position vector at the point s = t on the boundary, whilst r is the position 
vector inside the billiard B. G(r,r(t)) is the free particle Green function, namely 



G(r,r') = -^ 1) (A:|r-r'|), (15) 

where H$ (x) is the zero order Hankel function of the first kind. It is important to 
know that knowing Uj(s) for a certain eigenfunction j, we can immediately calculate 
the wave function ipj(r) within the interior of the billiard domain by 

^.( r ) = -£dtuj(t) G(r,r(t)). (16) 

Thus, in certain analogy to the classical mechanics, the quantum mechanics is completely 
described by the boundary functions Uj(s). Finally, there is the important identity [55] 

^jdt n(t) • r(t) Uj(t) 2 = k). (17) 

The quantum analogy of the classical phase space is the space of Wigner functions 
[55] or of other phase space representations of quantum states in general. The Wigner 
functions are real but not positive definite functions, and exhibit lots of oscillations 
around the zero level also in the regime where the quantum probability density is low 
and thus such structures often obscure the main physical aspects of the phenomena. 
Nevertheless, they do uniformly condense on the classical invariant objects, according 
to the mentioned PUSC [3] in the introduction. There are different ways of defining 
positive definite phase space functions, but Husimi functions [57] are perhaps the best 
way to do it. They are in fact Gaussian smoothed Wigner functions. In general, we 
define them by the projection of the wave function onto a coherent state. One of the 
possible formulations can be found in [M], whose definitions and the notation we shall 
use in what follows. The idea and the approach (in slightly different form) goes back to 
the works [58], [59] and [60]. The most important idea is to define the one-dimensional 
coherent states onto which we project the boundary functions Uj(s). For this reason, and 
due to the analogy with the classical dynamics and its Poincare surface of section on the 
cylinder (s,p), the underlying Husimi functions are called Poincare Husimi functions 

The key idea is to introduce one- dimensional coherent state as a function of the 
coordinate s on the boundary dB, localized at (q,p) G [0, C] x R, which is properly 
periodized, as introduced by Tualle and Voros [59J, but here following the notation from 
[541 we define 



exp{i kp(s — q + mC)} exp ^— — (s — q + mCf^j . (18) 



c (q,p)M S ) = 

The periodicity in s with the period C is now obvious. Here we have dropped all 
normalization factors, because in the end we shall normalize the Poincare Husimi 
functions anyway. Then, using this, the Poincare Husimi function associated with the 
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j-th eigenstate represented by the boundary function Uj(s) with the eigenvalue k = kj, 
is 



which is positive definite by construction. In the semiclassical limit j — > oo, and 
kj — > oo, we shall observe that the Poincare Husimi function is concentrated on the 
classical invariant regions, which can be an invariant torus, a chaotic component, or the 
entire Poincare surface of section (s,p) if the motion is ergodic. This is a consequence of 
PUSC, bearing in mind that the Husimi function is just a Gaussian smoothed Wigner 
function, where in the semiclassical limit the width of the smoothing Gaussian becomes 
less and less important. Therefore, we expect that in the semiclassical limit the Poincare 
Husimi functions will directly correspond to either the classical regular regions or to the 
classical chaotic regions, with the exceptions having measure zero. We can then use the 
Poincare Husimi functions to classify and thus also to separate the regular and chaotic 
eigenstates, and thereby also separate the regular and the chaotic spectral subsequences 
of the energy eigenvalues Ej = k 2 -. This is what we do in the next section |3j 

3. Separating the regular and chaotic eigenstates and subspectra 

We consider the billiard ( fTTT) with A = 0.15. As mentioned, the value of the classical 
parameter is p\ = 0.175. The quantum eigenstates were calculated using the method of 
Vergini and Saraceno [61] with great accuracy, for all eigenstates (587654) within the 
interval k G [2000, 2500]. The number of eigenstates below k = 2000 is estimated by the 
Weyl rule as about 1.000.000. Then, when calculating the Poincare Husimi functions, the 
momentum p is rescaled by the eigenvalue kj , such that p = 1 corresponds to the original 
p = kj. For each eigenstate the Poincare Husimi function was calculated as follows. We 
have set up a grid of 400 x 400 cells on the 1/4 of the surface of section (q,p), thus 
reduced due to the symmetries (reflection symmetry and time reversal symmetry). The 
grid points are defined as (qi,Pj) = (Aq/2 + iAq, Ap/2 + jAp), where Aq = £/800 and 
Ap = 1/400. The grid covers the 1/4 of the entire surface of section, namely q G [0, C/2] 
and p G [0, 1]. The grid points are positioned at the centers of the square cells of the 
area AqAp. The integration method used to evaluate (fl~9l) is a simple trapeze rule with 
the step ds ex A/20, where A = 2n/kj is the de Broglie wavelength. The important point 
is now that the values of the Husimi functions on the grid are normalized in such a way 
that their sum is equal to one. 

Some examples of the Poincare Husimi functions are shown in figure [TJ 
Now the classification of eigenstates can be performed by their projection onto the 
classical surface of section. As we are very deep in the semiclassical regime we do expect 
with probability one that either an eigenstate is regular or chaotic, with exceptions 
having measure zero, ideally. To automate this task we have ascribed to each point on 
the grid a number A^j whose value is either +1 if the grid point lies within the classical 
chaotic region or —1 if it belongs to a classical regular region. The Poincare Husimi 




(19) 
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Figure 1. Examples of chaotic (left) and regular (right) states in the Poincare-Husimi 
representation, kj (M) from top down are: chaotic: kj (M) = 2000.0021815 (0.978), 
2000.0181794 (0.981), 2000.0000068 (0.989), 2000.0258600 (0.965); regular: k 3 (M) = 
2000.0081402 (-0.987), 2000.0777155 ( -0.821), 2000.0786759 ( -0.528), 2000.0112417 ( 
-0.829). The gray background is the classically chaotic invariant component. We show 
only one quarter of the surface of section (s,p) £ [0,£/2] x [0, 1], because due to the 
reflection symmetry and time- reversal symmetry the four quadrants are equivalent. 



function H(q,p) (jl~9l) (normalized) was calculated on the grid points and the overlap 
index M was calculated according to the definition 

M y^ll,., A,. r (20) 

In practice, M is not exactly +1 or —1, but can have a value in between. The reasons 
are two, first the finite discretization of the phase space (the finite size grid), and second, 
the finite wavelength (not sufficiently small effective Planck constant, for which we can 
take just 1/kj). If so, the question is, where to cut the distribution of the M- values, 
at the threshold value M t , such that all states with M < M t are declared regular and 
those with M > M t chaotic. 

There are two natural criteria: (I) The classical criterion: the threshold value M t is 
chosen such that we have exactly p\ fraction of regular levels and P2 = 1 — p\ of chaotic 
levels. (II) The quantum criterion: we choose M t such that we get the best possible 
agreement of the chaotic level spacing distribution with the Brody distribution ([6]), which 
is expected to capture the dynamical localization effects of the chaotic eigenstates. 
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Figure 2. (a; left) The level spacing distribution for the entire spectrum after 
unfolding for N = 587653 spacings, with kj g [2000,2500], in excellent agreement 
with the BRB distribution with the classical p\ = 0.175 and (3 — 0.45. In the U- 
function plot (b; right), we show U(data) — U(BRB) as a function of S, and it is 
clearly seen that the BRB distribution with the quantally determined p\ = 0.19 and 
(3 = 0.47 is even better fit to the data (dashed, denoted by BRBq). The belt around 
the data curve indicates the expected statistical ± one-sigma error. 



4. Results 

To begin with we first look at the total energy spectrum Ej = for kj G [2000, 2500]. 
The spectral unfolding was done using the Weyl formula with the perimeter corrections. 
In figure [2^, we show the histogram of N = 587653 level spacings together with the 
best fitting BRB distribution ([2]), with E r being Poissonian and E c being the Brody gap 
probability ([8]), derived from ([6]), with the classical p\ = 0.175, and the Brody parameter 
(3 = 0.45. For the reference the BR distribution ([5]), the Poisson and the GOE level 
spacing distributions ([3]) are shown. In figure [2b we show the [/-function representation 
of the level spacing distribution, as introduced by Prosen and Robnik [23] and defined 
in the Appendix A. It was the [/-function which was used in finding the best fitting 
distributions. Whilst in the first case the agreement is perfect, in the [/-plot we see that 
the quantally adjusted parameter p\ = 0.19 instead of its classical value 0.175 leads to 
even better agreement. In this case = 0.47. Please observe that the deviations here 
are already extremely small, so the significance of the best fit is of extreme importance. 

Let us now separate the regular and chaotic eigenstates and the corresponding 
eigenvalues, after unfolding, according to the method described in section [3j using the 
classical criterion (I). The corresponding threshold value of the index M is found to be 
M t = 0.431. The level spacing distributions are shown in figure [3j As we see, we have 
perfect Brody distribution with /3 = 0.444 for the chaotic levels and almost pure Poisson 
for the regular levels. 

In order to make our analysis deeper and more refined we show first the histogram 



Dynamical localization in chaotic systems 12 



Chaotic levels Regular levels 




01234 01234 
S S 



Figure 3. Separation of levels using the classical criterion M t = 0.431. (a; left) 
The level spacing distribution for the chaotic subspectrum after unfolding, in perfect 
agreement with the Brody distribution (3 = 0.444. (b; right) The level spacing 
distribution for the regular part of the spectrum, after unfolding, in excellent agreement 
with Poisson. 



of the M- values in figure UK, where we see the two different threshold values M t , namely 
the classical one at M t = 0.431, corresponding to p\ = 0.175, and the quantum one at 
M t = 0.75 in which case we get p\ = 0.223. In figure Hb we analyze the goodness of the 
two fits, Brody and Izrailev. First we define the fit deviation measure as 

POO 

R = / P(S) 2 ((P(S) - P F (S)) 2 dS, (21) 
Jo 

where P(S) describes the data, whilst Pf(S) denotes the theoretical distribution fitting 
the data, namely F stands for Brody or Izrailev. At each chosen threshold value M t 
we consider the set of chaotic levels for which by definition M > M t . By performing 
the best fit at such M t (= threshold M) both for Brody and Izrailev, we calculate the 
fit deviation measure ( 121]) and plot its decadic logarithm as a function of M t in figure 
Hb- We see that at low M t Izrailev is somewhat better than Brody, but this is the 
unphysical domain of much too small M t . Near the classical threshold M t = 0.431 they 
become comparably good, but at increasing the M t Brody exhibits a very sharp and 
narrow minimum, whilst Izrailev curve increases. This sharp minimum of R for Brody 
is at the value of M which by definition we called the quantum threshold M t = 0.75. 
In addition, we show the R quantity ( 1211) for the Poisson distributioin vs. M t , where 
we see also a deep sharp minimum at M t « 0.8, thus almost at the quantum threshold 
M t = 0.75. The conclusion of this analysis is that by varying the threshold value M t 
there is a point M t = 0.75 at which the Brody distribution for all chaotic levels with 
M > M t is globally the best and also better than Izrailev fit at any M t . For logical 
consistency, it is important that at (almost) the same value of M t the fit of the Poisson 
distribution for the regular levels with M < M t is globally the best, as is evident from 
the R plot in figure Hb- 
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Figure 4. (a; left) Distribution of the index M with the locations of the threshold 
values of M, the classical M t = 0.431 and the quantum one M t = 0.75. In the first 
case classical p\ = 0.175, whilst in the second case the quantal p% — 0.223. (b; right) 
The logarithm of the fit deviation measure log 10 R, defined in (|2~Tj) . for the Brody and 
Izrailev distributions for chaotic levels M > M t , versus M t , and the same for Poisson 
distribution for the regular levels with M < M t . 




Figure 5. The [/-function plots as differences U(data) — U (ideal) for the regular and 
chaotic levels, for both criteria, the classical one and the quantum one. The belts 
around the data lines indicate the expected statistical ± one-sigma errors. 



We show the ^/-function plots for the regular levels and the chaotic levels in figure 
|5l using the Brody distribution, for the classical and quantum criterions, corresponding 
to the two different values of M t explained above. In both cases we plot the difference 
U(data) — U (ideal), so that in case of perfect agreement the line would coincide with the 
abscissa. We clearly see that the quantum criterion yields noticable better agreement 
than the classical criterion. Again, it must be emphasized that the agreement is 
extremely good, and the deviations of U(data) from U (ideal) are very small numbers. 



Dynamical localization in chaotic systems 14 



Best fitting Brody Best fitting Izrailev 




0.01 0.2 1 2 3 0.01 0.2 1 2 3 

S S 



Figure 6. We show the [/-function plot for the chaotic levels, clearly showing that 
Brody distribution (dashed) is much better than Izrailev distribution. On the left we 
show the best fitting Izrailev distribution at the point where Brody is globally the 
best fit using the quantum threshold M t = 0.75. On the right we show the Izrailev fit 
at the point M t — —0.5 where it is globally the best, and the best Brody fit at the 
same M t = —0.5, showing that also there the Brody fit is better than Izrailev. This 
is impressive because the effects are small, and the statistical significance very high. 
The belt around the data line indicates the expected statistical ± one-sigma error. 



Finally, we should comment on the relevance of the Brody distribution. The Brody 
distribution (12} 03J (JHJ) still has no theoretical foundation, but it definitely captures 
correctly the effects of dynamical (Chirikov) localization of chaotic eigenstates. This 
has been recently confirmed by Manos and Robnik [39J in case of the kicked rotator, 
namely for the quasienergies, and clearly is demonstrated in the present work for the 
autonomous Hamilton system, exemplified by the 2D billiard that we have chosen for 
this analysis. It remains as an open theoretical problem to derive the Brody distribution 
in this context. 

Since the Brody distribution is not known theoretically to be the preferred and/or 
"the right one" , we have considered again also the Izrailev distribution [35j |38j |3l] (JTUjl , 
studied very recently also in [39]. Entirely in line with the findings in [39] we found that 
Brody is much better model of the level spacing distribution than the Izrailev's one. 
This is clearly demonstrated in figure |6j 

5. Conclusions 

We have used a billiard system of the mixed type ( fTTT) . with A = 0.15, as introduced 
in [l5]-[l6], and have shown that using the Poincare Husimi functions we can separate 
the regular and chaotic eigenstates. The successful separation of course also entirely 
confirms the Berry-Robnik picture [5] of separating the regular and chaotic levels in 
the semiclassical limit, where the tunneling effects can be neglected. With great and 
unprecedented statistical significance we have shown that the chaotic levels exhibit 
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Brody level spacing distribution, whilst the regular levels obey Poissonian statistics. 
This analysis not only confirms the Berry-Robnik picture [5] of conceptually separating 
the regular and chaotic levels, based on the PUSC and embodied in formula ((H), but 
also demonstrates that the dynamical localization effects of the chaotic eigenstates are 
very well captured by the Brody distribution, in analogy with the same finding in the 
Floquet systems, in particular the kicked rotator [3UI39], where the quasienergy spectra 
are analyzed. Whilst in the kicked rotator the relationship between the localization 
measure of the eigenstates and the spectral level repulsion (Brody) parameter (3 exists, 
as proposed by Izrailev, and confirmed by Manos and Robnik, in the time independent 
Hamilton systems, like the one discussed in the present work, such relationship is 
lacking and is open for the future work. It involves great numerical efforts. The 
theoretical derivation of the Brody level spacing distribution for the dynamically 
localized eigenstates is also an open problem for the future. The billiard systems are not 
just nice theoretical toy models, but are suitable also for the experimental applications, 
like in quantum dots, and microwave cavities introduced and studied extensively over 
decades by H.-J. Stockmann []]. We also propose to study from the present point of 
view the hydrogen atom in strong magnetic field as an example of classical and quantum 
chaos par excellence, as introduced in [62J, [631 EU [651 E§] , although the technical efforts 
to obtain large stretches of high-lying eigenstates and the corresponding energy levels 
are much bigger than in billiard systems, where we have a great number of different 
elegant numerical techniques [67] , all of them used in our recent work [23] . 
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Appendix A: The U-function representation of the level spacing 
distribution 

First we estimate the expected fluctuation (error) of the cumulative (integrated) level 
spacing distribution W(S), which contains N s objects. At a certain S we have the 
probability W that a level is in the interval [0, W] and 1 — W that it is in the interval 
[W, 1]. Assuming binomial probability distribution P(k) of having k levels in the first 
and N q — k levels in the second interval we have 



W) 



N s -k 



(22) 



Then the average values are equal to 

< k >= N S W, < k 2 >= 



N S W + N S (N S 



l)W\ 



(23) 



and the variance 



V(k) =<k 2 > - <k> 2 



N S W(1 — W). 



(24) 
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But the probability W is estimated in the mean as k/N s . Its variance is 

and therefore the estimated error of W (standard deviation, the square root of the 
variance) is given by 



5W = VV{W) = \I W{ \ T W) - (26) 



N s 

Transforming now from W(S) to 

U(S) = - arccos y/l-W(S), (27) 

7T 

we show in a straightforward manner that 

SU = — ^= (28) 



and is indeed independent of S. From the (choice of the constant pre-factor in the) 
definition (|27l) one sees that both U(S) and W(S) go from to 1 as S goes from to 
infinity. 
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